In computer graphics, a computer is used to generate digital data that represents the projection of surfaces of objects in, for example, a three-dimensional scene, illuminated by one or more light sources, onto a two-dimensional image plane, to simulate the viewing of the scene by the eye of, for example, a person or the recording of the scene by, for example, a camera. The camera may include a lens for projecting the image of the scene onto the image plane, or it may comprise a pinhole camera in which case no lens is used. The two-dimensional image is in the form of an array of picture elements (which are variable termed “pixels” or “Pels”), and the digital data generated for each pixel represents the color and luminance of the scene as projected onto the image plane at the point of the respective pixel in the image plane. The surfaces of the objects may have any of a number of characteristics, including shape, color, specularity, texture, and so forth, which are preferably rendered in the image as closely as possible, to provide a realistic-looking image.
Generally, the contributions of the light reflected from the various points in the scene to the pixel value representing the color and intensity of a particular pixel are expressed in the form of the one or more integrals of relatively complicated functions. Since the integrals used in computer graphics generally will not have a closed-form solution, numerical methods must be used to evaluate them and thereby generate the pixel value. Typically, a conventional “Monte Carlo” method has been used in computer graphics to numerically evaluate the integrals. Generally, in the Monte Carlo method, to evaluate an integral                               〈          f          〉                =                              ∫            0            1                    ⁢                                    f              ⁡                              (                x                )                                      ⁢                                                   ⁢                          ⅆ              x                                                          (        1        )            where f(x) is a real function on the “s”-dimensional unit cube [0,1)s (that is, an s-dimensional cube that includes “zero,” and excludes “one”), first a number “N” statistically-independent randomly-positioned points xi , i=1, . . . , N, are generated over the integration domain. The random points xi are used as sample points for which sample values f(xi) are generated for the function f(x), and an estimate {overscore (ƒ)} for the integral is generated as                                           〈            f            〉                    ≈                      f            _                          =                              1            N                    ⁢                                    ∑                              i                =                1                            N                        ⁢                                                   ⁢                          f              ⁡                              (                                  x                  i                                )                                                                        (        2        )            As the number of random points used in generating the sample points f(xi) increases, the value of the estimate {overscore (ƒ)} will converge toward the actual value of the integral <ƒ>. Generally, the distribution of estimate values that will be generated for various values of “N,” that is, for various numbers of sample points, of being normal distributed around the actual value with a standard deviation which can be estimated by                     σ        =                                            1                              N                -                1                                      ⁢                          (                                                                    f                    2                                    _                                -                                                      f                    _                                    2                                            )                                                          (        3        )            if the points xi used to generate the sample values f(xi) are statistically independent, that is, if the points xi are truly positioned at random in the integration domain.
Generally, it has been believed that random methodologies like the Monte Carlo method are necessary to ensure that undesirable artifacts, such as Moiré patterns and aliasing and the like, which are not in the scene, will not be generated in the generated image. However, several problems arise from use of the Monte Carlo method in computer graphics. First, since the sample points xi used in the Monte Carlo method are randomly distributed, they may clump in various regions over the domain over which the integral is to be evaluated. Accordingly, depending on the set of points that are generated, in the Monte Carlo method for significant portions of the domain there may be no sample points xi for which sample values f(xi) are generated. In that case, the error can become quite large. In the context of generating a pixel value in computer graphics, the pixel value that is actually generated using the Monte Carlo method may not contain some elements which might otherwise be contained if the sample points xi were guaranteed to be more evenly distributed over the domain. This problem can be alleviated somewhat by dividing the domain into a plurality of sub-domains, but it is generally difficult to determine a priori the number of sub-domains into which the domain should be divided, and, in addition, in a multi-dimensional integration domain, the partitioning of the integration domain into sub-domains, which are preferably of equal size, can be quite complicated.
In addition, since the method makes use of random numbers, the error |{overscore (ƒ)}−<ƒ>| (where |x| represents the absolute value of the value “x”) between the estimate value {overscore (ƒ)} and actual value <ƒ> is probabilistic, and, since the error values for various large values of “N” are close to normal distribution around the actual value <ƒ>, only sixty-eight percent of the estimate values {overscore (ƒ)} that might be generated are guaranteed to lie within one standard deviation of the actual value <ƒ>.
Furthermore, as is clear from equation (3), the standard deviation decreases with increasing numbers “N” of sample points, proportional to the reciprocal of square root of “N”(that is, √{square root over (1/N)}).
Thus, if it is desired to reduce the statistical error by a factor of two, it will be necessary to increase the number of sample points N by a factor of four, which, in turn, increases the computational load that is required to generate the pixel values, for each of the numerous pixels in the image.
Additionally, since the Monte Carlo method requires random numbers to define the coordinates of respective sample points xi in the integration domain, an efficient mechanism for generating random numbers is needed. Generally, digital computers are provided with so-called “random number” generators, which are computer programs which can be processed to generate a set of numbers that are approximately random. Since the random number generators use deterministic techniques, the numbers that are generated are not truly random. However, the property that subsequent random numbers from a random number generator are statistically independent should be maintained by deterministic implementations of pseudo-random numbers on a computer.
The Grabenstein application describes a computer graphics system and method for generating pixel values for pixels in an image using a strictly deterministic methodology for generating sample points, which avoids the above-described problems with the Monte Carlo method. The strictly deterministic methodology described in the Grabenstein application provides a low-discrepancy sample point sequence that ensures, a priori, that the sample points are generally more evenly distributed throughout the region over which the respective integrals are being evaluated. In one embodiment described in the Grabenstein application, the sample points that are used are based on a so-called Halton sequence. See, for example, J. H. Halton, Numerische Mathematik, Vol. 2, pp. 84-90 (1960) and W. H. Press, et al., Numerical Recipes in Fortran (2d Edition) page 300 (Cambridge University Press, 1992). In a Halton sequence generated for number base “b,” where base “b” is a selected prime number, the “i-th” value of the sequence, represented by Hbi is generated by use of a “radical inverse” function Φb that is generally defined as                                                         Φ              b                        ⁢                          :                        ⁢                                                   ⁢                          N              0                                →          I                ⁢                                  ⁢                  i          =                                                    ∑                                  j                  =                  0                                ∞                            ⁢                                                           ⁢                                                                    a                    j                                    ⁡                                      (                    i                    )                                                  ⁢                                  b                  j                                                      ↦                                          ∑                                  j                  =                  0                                ∞                            ⁢                                                           ⁢                                                                    a                    j                                    ⁡                                      (                    i                    )                                                  ⁢                                  b                                                            -                      j                                        -                    1                                                                                                          (        4        )            where (αj)j=0∞ is the representation of “i” in integer base “b.” Generally, a radical inverse of a value “k” is generated by                (1) writing the value “i” as a numerical representation of the value in the selected base “b,” thereby to provide a representation for the value as DM DM-1 . . . D2D1, where Dm(m=1,2, . . . , M) are the digits of the representation,        (2) putting a radix point (corresponding to a decimal point for numbers written in base ten) at the least significant end of the representation DMDM-1 . . . D2D1 written in step (1) above, and        (3) reflecting the digits around the radix point to provide 0.D1D2 . . . DM-1DM, which corresponds to Hbi.It will be appreciated that, regardless of the base “b” selected for the representation, for any series of values, one, two, . . . “k,” written in base “b,” the least significant digits of the representation will change at a faster rate than the most significant digits. As a result, in the Halton sequence Hb1, Hb2, . . . Hbi, the most significant digits will change at the faster rate, so that the early values in the sequence will be generally widely distributed over the interval from zero to one, and later values in the sequence will fill in interstices among the earlier values in the sequence. Unlike the random or pseudo-random numbers used in the Monte Carlo method as described above, the values of the Halton sequence are not statistically independent; on the contrary, the values of the Halton sequence are strictly deterministic, “maximally avoiding” each other over the interval, and so they will not clump, whereas the random or pseudo-random numbers used in the Monte Carlo method may clump.        
It will be appreciated that the Halton sequence as described above provides a sequence of values over the interval from zero to one, inclusive along a single dimension. A multi-dimensional Halton sequence can be generated in a similar manner, but using a different base for each dimension.
A generalized Halton sequence, of which the Halton sequence described above is a special case, is generated as follows. For each starting point along the numerical interval from zero to one, inclusive, a different Halton sequence is generated. Defining the pseudo-sum x⊕py for any x and y over the interval from zero to one, inclusive, for any integer “p” having a value greater than two, the pseudo-sum is formed by adding the digits representing “x” and “y” in reverse order, from the most-significant digit to the least-significant digit, and for each addition also adding in the carry generated from the sum of next more significant digits. Thus, if “x” in base “b” is represented by 0.X1X2 . . . XM-1XM, where each “Xm” is a digit in base “b,” and if “y” in base “b” is represented by 0.Y1Y2 . . . YN-1YN, where each “Yn” is a digit in base “b” (and where “M,” the number of digits in the representation of “x” in base “b”, and “N,” the number of digits in the representation of “y” in base “b”, may differ), then the pseudo-sum “z” is represented by 0.Z1Z2 . . . ZL-1ZL, where each “Z1” is a digit in base “b” given by Zl=(Xl+Yl+Cl) mod b, where “mod” represents the modulo function, and       C    l    =      {                            1                                                                    for                ⁢                                                                   ⁢                                  X                                      l                    -                    1                                                              +                              Y                                  l                  -                  1                                            +                              Z                                  l                  -                  1                                                      ≥            b                                                0                          otherwise                    is a carry value from the “1-1st” digit position, with Cl being set to zero.
Using the pseudo-sum function as described above, the generalized Halton sequence that is used in the system described in the Grabenstein application is generated as follows. If “b” is an integer, and x0 is an arbitrary value on the interval from zero to one, inclusive, then the “p”-adic von Neumann-Kakutani transformation Tb(x) is given by                                                         T              p                        ⁡                          (              x              )                                :=                      x            ⁢                          ⊕              p                        ⁢                          1              b                                      ,                            (        5        )            and the generalized Halton sequence x0, x1, x2, . . . is defined recursively asxn+1=Tb(xn)  (6)From equations (5) and (6), it is clear that, for any value for “b,” the generalized Halton sequence can provide that a different sequence will be generated for each starting value of “x,” that is, for each x0. It will be appreciated that the Halton sequence Hbi as described above is a special case of the generalized Halton sequence (equations (5) and (6)) for x0=0.
The Keller application also describes a computer graphics system and method for generating pixel values for pixels in an image using a strictly deterministic methodology, in this case using s-dimensional Hammersley point sets. S-dimensional Hammersley point sets are defined as                                                                                           U                                      N                    ,                    s                                    Hammersley                                ⁢                                  :                                ⁢                                                                   ⁢                                  {                                      0                    ,                    …                    ⁢                                                                                   ,                                          N                      -                      1                                                        }                                            →                            ⁢                              I                s                                                                                                                          i                  ↦                                    ⁢                                      x                    i                                                  :=                                  (                                                            i                      N                                        ,                                                                  Φ                                                  b                          1                                                                    ⁡                                              (                        i                        )                                                              ,                    …                    ⁢                                                                                   ,                                                                  Φ                                                  b                                                      s                            -                            1                                                                                              ⁡                                              (                        i                        )                                                                              )                                            ,                                                          (        7        )            where Is is the s-dimensional unit cube [0,1), the number of points “N” in the set is fixed and Φbx(i) is the radical inverse of “i” in base bx, for bases b1, . . . ,bs-1, as sample points. The bases do not need to be prime, but they are preferably relatively prime to provide a relatively uniform distribution.
In addition, the Keller application also describes use of a scrambled Hammersley point set to reduce or eliminate a problem that can arise in connection with higher-dimensioned low-discrepancy sequences. Typically, the radical inverse function Φb has subsequences of b−1 equidistant values spaced by 1/b, which can result in correlation patterns. Although these correlation patterns are merely noticeable in the full s-dimensional space, they are undesirable since they are prone to aliasing. The Keller application describes a methodology that attenuates this effect by scrambling, which corresponds to application of a permutation to the digits of the b-ary representation used in the radical inversion. For a permutation a from a symmetric group Sb over integers 0, . . . ,b−1, the scrambled radical inverse is defined as                                                                                           Φ                  b                                ⁢                                  :                                ⁢                                                                   ⁢                                  N                  0                                ×                                  S                  b                                            →                            ⁢              I                                                                                                            (                                      i                    ,                    σ                                    )                                ↦                                ⁢                                                      ∑                                          j                      =                      0                                        ∞                                    ⁢                                                                           ⁢                                                            σ                      ⁡                                              (                                                                              a                            j                                                    ⁡                                                      (                            i                            )                                                                          )                                                              ⁢                                          b                                                                        -                          j                                                -                        1                                                                                            ⇔                i                            =                                                ∑                                      j                    =                    0                                    ∞                                ⁢                                                                   ⁢                                                                            a                      j                                        ⁡                                          (                      i                      )                                                        ⁢                                                            b                      j                                        .                                                                                                          (        8        )            If the permutation “σ” is the identity, the scrambled radical inverse corresponds to the unscrambled radical inverse. In one embodiment, the permutation σ is defined recursively as follows. Starting from the permutation σ2=(0,1) for base b=2, the sequence of permutations is defined as follows:                (i) if the base “b” is even, the permutation σb is generated by first taking the values of   2  ⁢      σ          b      2      and appending the values of             2      ⁢              σ                  b          2                      +    1    ,and        (ii) if the base “b” is odd, the permutation σb is generated by taking the values of σb-1, incrementing each value that is greater than or equal to       b    -    1    2by one, and inserting the value b−1 in the middle.        
This recursive procedure results in the following permutations                σ2=(0,1)        σ3=(0,1,2)        σ4=(0,2,1,3)        σ5=(0,3,2,1,4)        σ6=(0,2,4,1,3,5)        σ7=(0,2,5,3,1,4,6)        σ8=(0,4,2,6,1,5,3,7) . . .Accordingly, given the radical inverse for a value “i” in base “b,” if the “k-th” digit of the representation for the radical inverse has the value “j,” the “k-th” digit of the scrambled radical inverse has the value corresponding to the value of the “j-th” digit in the permutation σb above. Use of the scrambled radical inverse, which may be applied to Halton and other low-discrepancy sequences as well as the Hammersley point set, reduces or eliminates the correlation patterns described above.        
The use of a strictly deterministic low-discrepancy sequence can provide a number of advantages over the random or pseudo-random numbers that are used in connection with the Monte Carlo technique. Unlike the random numbers used in connection with the Monte Carlo technique, the low discrepancy sequences ensure that the sample points are more evenly distributed over a respective region or time interval, thereby reducing error in the image which can result from clumping of such sample points which can occur in the Monte Carlo technique. That can facilitate the generation of images of improved quality when using the same number of sample points at the same computational cost as in the Monte Carlo technique. However, problems can arise if it is desired to generate an image using a plurality of processors operating in parallel, since it will generally be desirable to ensure that each of the processors makes use of a low discrepancy sequence in its processing operations.